#import <sys/message.h>
#import <mach_error.h>
#import <soundkit/Sound.h>
#import <string.h>
#import "Pvc_Object.h"

void bitreverse( float *buf,int N );

@implementation Pvc_Object : Object
{

/*

  int		N,
		N2,
		Nw = 0,
		Nw2, 
		decim = 0,
		interp = 0,
		in,
		on,
		valid = -1,
		eof = 0;
  float 	freqmlt = 0.,
		srate = 44100.,
		pi,
		twopi,
		synt = 0.,
		*Hwin,
		*Wanal,
		*Wsyn,
		*input,
		*buf,
		*channel,
		*output;
  BOOL		obank = 0,
		aflag = 0,
		sflag = 0;
  NXZone	*ozone;

*/

}

+ create {
    id 		newInstance;
    newInstance = [ self new ];
    return newInstance;
}

- init {
    [super init];
    [self setDefaults];
    inputFilename = malloc(1024);
    status = IDLE;
    return self;
}

- setInputFilename: (char *) arg
{
    int err;
    if (inputFilename && !(strcmp(arg,inputFilename)))
    	return self;
    if (insfh != NULL)
	SNDFree(insfh);    	
    if (strlen(arg) > 1023) 
    	return nil;
    strcpy(inputFilename,arg);
    err = SNDReadSoundfile(inputFilename, &insfh);
    if (err != SND_ERR_NONE)
	return nil;
    return self;
}

- setOutputFilename: (char *) arg
{
    outputFilename = arg;
    return self;
}

static BOOL isPowerOfTwo(int n)
    /* Query whether n is a pure power of 2 */
{
    while (n > 1) {
	if (n % 2) break;
	n >>= 1;
    }
    return (n == 1);
}

- checkArgs:(char *)msg 
 /* If failure, sets msg to error message and returns nil.
    Msg should be at least 1024 long */
{

    if (Nw == 0)
	Nw = N;

    if (interp == 0)
      interp = decim;

    obank = freqmlt != 0.;
    N2 = N>>1;
    Nw2 = Nw>>1;
    if (interp > Nw) {
	sprintf(msg,"I must be less than or equal to window size.");
	return nil;
    }
    if (decim > Nw) {
	sprintf(msg,"D must be less than or equal to window size.");
	return nil;
    }
    if (decim > Nw) {
	sprintf(msg,"D must be less than or equal to window size.");
	return nil;
    }
    if (!isPowerOfTwo(N)) {
	sprintf(msg,"FFT size (N) must be a power of 2.");
    }
    return self;
}

- allocateDataspace {

    int		width,
		size;

    if ( ozone != NULL )
      NXDestroyZone(ozone);

    ozone = NXCreateZone( 10 * Nw * sizeof(float), vm_page_size, 0 );

/* analysis window */
    Wanal = (float *) zspace( ozone, Nw, sizeof(float) );

 /* synthesis window */
    Wsyn = (float *) zspace( ozone, Nw, sizeof(float) );

/* input buffer */
    input = (float *) zspace( ozone, Nw, sizeof(float) );

/* hamming window */
    Hwin = (float *) zspace( ozone, Nw, sizeof(float) );

/* FFT buffer */
    buf = (float *) zspace( ozone, N, sizeof(float) );

/* analysis channels */
    channel = (float *) zspace( ozone, N+2, sizeof(float) );

/* output buffer */
    output = (float *) zspace( ozone, Nw, sizeof(float) );

    SNDGetDataPointer( insfh, (char **)&inData, &insize, &width );
    if (outsfh != NULL)
    	SNDFree(outsfh);
    SNDAlloc( &outsfh, 
	     ( (int) ((N+insize) * width * ( (float) interp / decim ))),
    		SND_FORMAT_LINEAR_16, insfh->samplingRate,
	        insfh->channelCount, 0);
    SNDGetDataPointer(outsfh, (char **)&outData, &size, &width );

    srate = insfh->samplingRate;

    return self;
}

- (SNDSoundStruct *)outsfh
{
  return outsfh;
}  

- (SNDSoundStruct *)insfh
{
  return insfh;
}  

- setDefaults {

  N = 1024;
  Nw = 0;
  valid = -1;
  decim = 128; 
  interp = 0;
  eof = 0;
  freqmlt = 0.;
  srate = 44100.;
  pi = 4 * atan(1.);
  twopi = 8 * atan(1.);
  synt = 0.;
  obank = 0;
  inPoint = 0;
  outPoint = 0;

  return self;
}

- (float)howFar {
  return ((float)outPoint) / ((outsfh->dataSize)/sizeof(short));
}

- N: (int) aN {
    N = aN;
    return self;
}

- Nw: (int) aNw {
    Nw = aNw;
    return self;
}

- decim: (int) aDecim {
    decim = aDecim;
    return self;
}

- interp: (int) aInterp {
    interp = aInterp;
    return self;
}

- freqmlt: (float) aFreqmlt {
    freqmlt = aFreqmlt;
    return self;
}

- srate: (float) aSrate {
    srate = aSrate;
    return self;
}

- synt: (float) aSynt {
    synt = aSynt;
    return self;
}

- aflag: (BOOL) aAflag {
    aflag = aAflag;
    return self;
}

- sflag: (BOOL) aSflag {
    sflag = aSflag;
    return self;
}

- writeOutputSound {

  int err;
  err = SNDWriteSoundfile( outputFilename, outsfh );
  if (err != SND_ERR_NONE)
  	return nil;
  return self;
}

- kill
{
    [self resume];
    eof = 1;
    return self;
}

- resume
  /* Gets called from the other thread */
{
    kern_return_t ec;
    msg_header_t msg =    {0,                   /* msg_unused */
                           TRUE,                /* msg_simple */
                           sizeof(msg_header_t),/* msg_size */
                           MSG_TYPE_NORMAL,     /* msg_type */
                           0};                  /* Fills in remaining fields */
    if (eof != 2 && appToObjPort != PORT_NULL)
      return nil;
    msg.msg_local_port = PORT_NULL;
    msg.msg_remote_port = appToObjPort;
    ec = msg_send(&msg, SEND_TIMEOUT, 0);
    if (ec == SEND_TIMED_OUT)
      ; /* message queue is full, don't need to send another */
    return self;
}

- pause
{
    eof = 2;
    return self;
}

- (int)status
{
    return status;
}

- _waitForMessage
    /* This is the main loop for a separate-threaded performance. */
{
    struct {
        msg_header_t header;
        char data[MSG_SIZE_MAX];
    } msg;
    status = PAUSED;
    (void)port_allocate(task_self(), &appToObjPort);
    msg.header.msg_size = MSG_SIZE_MAX;
    msg.header.msg_local_port = appToObjPort;
    (void)msg_receive(&msg.header, MSG_OPTION_NONE,0);
    (void)port_deallocate(task_self(),appToObjPort);
    appToObjPort = PORT_NULL;
    status = RUNNING;
    return self;
  }

/* phase vocoder */

- runPvc
	/* Returns self if successfully writes file, else nil. */
{

    status = RUNNING;
/* initialize input and output time values (in samples) */
    in = -Nw;
    if ( decim )
	on = (in * interp) / decim;
    else
	on = in;

/* main loop--perform phase vocoder analysis-resynthesis */

    while ( eof != 1) {

        if (eof == 2)  {
	  [self _waitForMessage];
	  if (eof == 1)
	    break;
	}
        
/* increment times */
	in += decim;
	on += interp;

/* analysis: input D samples; window, fold and rotate input
   samples into FFT buf; take FFT; and convert to
   amplitude-frequency (phase vocoder) form */

	eof = [self shiftin];
	[self fold];
	[self rfft:YES];
	[self convert];

/* at this point channel[2*i] contains amplitude data and
   channel[2*i+1] contains frequency data (in Hz) for phase
   vocoder channels i = 0, 1, ... N/2; the center frequency
   associated with each channel is i*f, where f is the
   fundamental frequency of analysis R/N; any desired spectral
   modifications can be made at this point: pitch modifications
   are generally well suited to oscillator bank resynthesis,
   while time modifications are generally well (and more
   efficiently) suited to overlap-add resynthesis */

/* oscillator bank resynthesis */

	if ( obank ) {

	    [self oscbank]; 
/* oscbank( channel, N2, srate, Nw, interp, freqmlt, output ); */
	    [self shiftout];
	} 

/* overlap-add resynthesis */

	else {
	    [self unconvert];
	    [self rfft:NO];
	    [self overlapadd];
	    [self shiftout];
	}
    }
    if (![self writeOutputSound]) {
	status = IDLE;
	return nil;
    }
    status = IDLE;
    return self;
}


/* overlap-add FFT analysis/resynthesis without phase unwrapping */

- runFFT
{

/* argument checking */
    status = RUNNING;

/* initialize input and output time values (in samples) */
    in = -Nw;
    if ( decim )
	on = (in * interp) / decim;
    else
	on = in;

/* main loop--perform phase vocoder analysis-resynthesis */

    while ( !eof ) {

/* increment times */
	in += decim;
	on += interp;

/* analysis: input D samples; window, fold and rotate input
   samples into FFT buf; take FFT; and convert to
   amplitude-frequency (phase vocoder) form */

	eof = [self shiftin];
	[self fold];
	[self rfft:YES];
	[self convertFFT];
	
/* at this point channel[2*i] contains amplitude data and
   channel[2*i+1] contains frequency data (in Hz) for phase
   vocoder channels i = 0, 1, ... N/2; the center frequency
   associated with each channel is i*f, where f is the
   fundamental frequency of analysis R/N; any desired spectral
   modifications can be made at this point: pitch modifications
   are generally well suited to oscillator bank resynthesis,
   while time modifications are generally well (and more
   efficiently) suited to overlap-add resynthesis */

/* overlap-add resynthesis */

	[self unconvertFFT];
	[self rfft:NO];
	[self overlapadd];
	[self shiftout];
      }
    status = IDLE;
    return self;
}

/* shift next D samples into righthand end of array A of
   length N, padding with zeros after last sample (A is
   assumed to be initially 0); return 0 when more input
   remains, otherwise return 1 after N-2*D zeros have been
   padded onto end of input */

- (int) shiftin
{
  register int 	i;

  if ( valid < 0 )		/* first time only */
    valid = Nw;

  for ( i = 0 ; i < Nw - decim ; i++ )
    *(input+i) = *(input + i + decim);

  if ( valid == Nw ) {
    for ( i = (Nw - decim); i < Nw; i++) {
      if ( inPoint >= insize ) {
      	*(input+i) = (float) (*(inData + inPoint) / 32768.);
	valid = i - (Nw-decim);
	break;
      }
      *(input+i) = (float) (*(inData + inPoint) / 32768.);
      ++inPoint;
    }
  }
  if ( valid < Nw ) {		/* pad with zeros after EOF */
    for (i=valid; i < Nw; i++)
      *(input+i) = 0.;
    valid -= decim;
  }
  return (valid <= 0);
}


/* if output time n >= 0, output first I samples in
   array A of length N, then shift A left by I samples,
   padding with zeros after last sample */

- shiftout
{
  int	i;

  if ( on >= 0 ) {
    for ( i=0; i < interp; i++ )
      *(outData + outPoint + i) = (short) (32767. * *(output+i));
      
    outPoint += interp;
  }

  for ( i = 0; i < Nw - interp; i++ )
    *(output+i) = *(output + i + interp);

  for ( i = (Nw - interp); i < Nw; i++ )
    *(output+i) = 0.;

  return self;
}

/* S is a spectrum in rfft format, i.e., it contains N real values
   arranged as real followed by imaginary values, except for first
   two values, which are real parts of 0 and Nyquist frequencies;
   convert first changes these into N/2+1 PAIRS of magnitude and
   phase values to be stored in output array C; the phases are then
   unwrapped and successive phase differences are used to compute
   estimates of the instantaneous frequencies for each phase vocoder
   analysis channel; decimation rate D and sampling rate R are used
   to render these frequency values directly in Hz. */

- convert
{
  static int 	first = 1;
  static float 	*lastphase,
		fundamental,
		factor;
  float 	phase,
		phasediff;
  int 		real,
		imag,
		amp,
		freq;
  float 	a,
		b;
  int 		i;

/* first pass: allocate zeroed space for previous phase
   values for each channel and compute constants */

    if ( first ) {
      first = 0;
      lastphase = (float *) space( N2+1, sizeof(float) );
      bzero(lastphase,(N2+1)*sizeof(float));
      fundamental = srate / (float) (N2<<1);
      factor = srate / (decim * twopi);
    } 

/* unravel rfft-format spectrum: note that N2+1 pairs of
   values are produced */

    for ( i = 0; i <= N2; i++ ) {
      imag = freq = ( real = amp = i<<1 ) + 1;
      a = ( i == N2 ? *(buf+1) : *(buf+real) );
      b = ( i == 0 || i == N2 ? 0. : *(buf+imag) );

/* compute magnitude value from real and imaginary parts */

      *(channel+amp) = hypot( a, b );

/* compute phase value from real and imaginary parts and take
   difference between this and previous value for each channel */

      if ( *(channel+amp) == 0. )
	phasediff = 0.;
      else {
	phasediff = ( phase = -atan2( b, a ) ) - *(lastphase+i);
	*(lastphase+i) = phase;
	
/* unwrap phase differences */

	while ( phasediff > pi )
	  phasediff -= twopi;
	while ( phasediff < -pi )
	  phasediff += twopi;
      }

/* convert each phase difference to Hz */

      *(channel+freq) = (phasediff * factor) + (i * fundamental);
    }
    return self;
}

- convertFFT
{

  int 		amp,
		phase;
  float 	a,b;
  register int	i;


    for ( i = 0; i <= N2; i++ ) {
      phase = (amp = i<<1) + 1;
      a = ( i == N2 ? *(buf+1) : *(buf+amp) );
      b = ( i == 0 || i == N2 ? 0. : *(buf+phase) );

/* compute magnitude value from real and imaginary parts */

      *(channel+amp) = hypot( a, b );
      *(channel+phase) = -atan2( b, a );
    }
    return self;
}


/* unconvert essentially undoes what convert does, i.e., it
  turns N2+1 PAIRS of amplitude and frequency values in
  C into N2 PAIR of complex spectrum data (in rfft format)
  in output array S; sampling rate R and interpolation factor
  I are used to recompute phase values from frequencies */

- unconvert
{
  static int 	first = 1;
  static float 	*lastphase,
		fundamental,
		factor;
  int 		i,
		real,
		imag,
		amp,
		freq;
  float 	phase;

/* first pass: allocate memory and compute constants */

    if ( first ) {
	first = 0;
	lastphase = (float *) space( N2+1, sizeof(float) );
	fundamental = srate / (float) (N2<<1);
        factor = twopi * interp / srate;
    } 

/* subtract out frequencies associated with each channel,
   compute phases in terms of radians per I samples, and
   convert to complex form */

    for ( i = 0; i <= N2; i++ ) {
	imag = freq = ( real = amp = i<<1 ) + 1;

	if ( i == N2 )
	    real = 1;

	*(lastphase+i) += *(channel+freq) - i * fundamental;
	phase = *(lastphase+i) * factor;
	*(buf+real) = *(channel+amp) * cos( phase );

	if ( i != N2 )
	    *(buf+imag) = - *(channel+amp) * sin( phase );
    }
    return self;
}


- unconvertFFT
{

  int 		amp,
		phase;
  register int	i;


    for ( i = 0; i <= N2; i++ ) {
	phase = (amp = i<<1) + 1;

	if ( i == N2 )
	  *(buf+1) = *(channel+amp) * cos( *(channel+phase) );
	else
	  *(buf+amp) = *(channel+amp) * cos( *(channel+phase) );

	if ( i != N2 )
	    *(buf+phase) = - *(channel+amp) * sin( *(channel+phase) );
    }
    return self;
}


/* make balanced pair of analysis (A) and synthesis (S) windows;
   window lengths are Nw, FFT length is N, synthesis interpolation
   factor is I, and osc is true (1) if oscillator bank resynthesis 
   is specified */

- makewindows
{
  int i;
  float sum;

/* basic Hamming windows */

  for ( i = 0; i < Nw; i++ )
    *(Hwin+i) = *(Wanal+i) = *(Wsyn+i) = 0.54 - (0.46 * cos( 
		twopi * i / (Nw - 1) ));

/* when Nw > N, also apply interpolating (sinc) windows to
   ensure that window are 0 at increments of N (the FFT length)
   away from the center of the analysis window and of I away
   from the center of the synthesis window */

    if ( Nw > N ) {
     float x;

/* take care to create symmetrical windows */

	x = -(Nw - 1) / 2.;
	for ( i = 0; i < Nw; i++, x += 1. )
	    if ( x != 0. ) {
		*(Wanal+i) *= N * sin( (pi * x) / N ) / (pi*x);
		if ( interp )
		   *(Wsyn+i) *= interp * sin( (pi*x) / interp ) / (pi*x);
	    }
    }

/* normalize windows for unity gain across unmodified
   analysis-synthesis procedure */

    for ( sum = i = 0; i < Nw; i++ )
	sum += *(Wanal+i);

    for ( i = 0; i < Nw; i++ ) {
     float afac = 2./sum;
     float sfac = Nw > N ? 1./afac : afac;
	*(Wanal+i) *= afac;
	*(Wsyn+i) *= sfac;
    }

    if ( Nw <= N && interp ) {
	for ( sum = i = 0; i < Nw; i += interp )
	    sum += *(Wsyn+i) * *(Wsyn+i);
	for ( sum = 1./sum, i = 0; i < Nw; i++ )
	    *(Wsyn+i) *= sum;
    }
    return self;
}


/* input I is a folded spectrum of length N; output O and
   synthesis window W are of length Nw--overlap-add windowed,
   unrotated, unfolded input data into output O */

- overlapadd
{
  int 	i,
	n;

  n = on;
  while ( n < 0 )
    n += N;
  n %= N;

  for ( i = 0 ; i < Nw ; i++ ) {
    *(output+i) += *(buf+n) * *(Wsyn+i);
    if ( ++n == N )
      n = 0;
  }
  return self;
}


/* multiply current input I by window W (both of length Nw);
   using modulus arithmetic, fold and rotate windowed input
   into output array O of (FFT) length N according to current
   input time n */

- fold
{
 
  int 	i,
	n;

  n = in;

  for ( i = 0; i < N; i++ )
    *(buf+i) = 0.;

    while ( n < 0 )
	n += N;

    n %= N;
    for ( i = 0; i < Nw; i++ ) {
      *(buf+n) += *(input+i) * *(Wanal+i);
      if ( ++n == N )
	n = 0;
    }
    return self;
}

/* If forward is true, rfft replaces 2*N real data points in buf with
   N complex values representing the positive frequency half of their
   Fourier spectrum, with *(buf+1) replaced with the real part of the Nyquist
   frequency value.  If forward is false, rfft expects buf to contain a
   positive frequency spectrum arranged as before, and replaces it with
   2*N real values.  N MUST be a power of 2. */

- rfft: (BOOL) forward
{
  float 	c1,c2,
  		h1r,h1i,
		h2r,h2i,
		wr,wi,
		wpr,wpi,
  		temp,
		theta;
  float 	br,bi;
  int 		i,
		i1,i2,i3,i4,
		N2p1;

    theta = pi/N2;
    wr = 1.;
    wi = 0.;
    c1 = 0.5;
    if ( forward ) {
	c2 = -0.5;
	[self cfft:forward];
	br = *buf;
	bi = *(buf+1);
    }
    else {
	c2 = 0.5;
	theta = -theta;
	br = *(buf+1);
	bi = 0.;
	*(buf+1) = 0.;
    }
    wpr = -2. * pow( sin( 0.5*theta ), 2. );
    wpi = sin(theta);
    N2p1 = (N2 <<1) + 1;
    for ( i = 0; i <= N2 >> 1; i++ ) {
	i1 = i<<1;
	i2 = i1 + 1;
	i3 = N2p1 - i2;
	i4 = i3 + 1;
	if ( i == 0 ) {
	    h1r =  c1 * ( *(buf+i1) + br );
	    h1i =  c1 * ( *(buf+i2) - bi );
	    h2r = -c2 * ( *(buf+i2) + bi );
	    h2i =  c2 * ( *(buf+i1) - br );
	    *(buf+i1) =  h1r + (wr * h2r) - (wi * h2i);
	    *(buf+i2) =  h1i + (wr * h2i) + (wi * h2r);
	    br =  h1r - (wr * h2r) + (wi * h2i);
	    bi = -h1i + (wr * h2i) + (wi * h2r);
	}
	else {
	    h1r =  c1 * ( *(buf+i1) + *(buf+i3) );
	    h1i =  c1 * ( *(buf+i2) - *(buf+i4) );
	    h2r = -c2 * ( *(buf+i2) + *(buf+i4) );
	    h2i =  c2 * ( *(buf+i1) - *(buf+i3) );
	    *(buf+i1) =  h1r + wr*h2r - wi*h2i;
	    *(buf+i2) =  h1i + wr*h2i + wi*h2r;
	    *(buf+i3) =  h1r - wr*h2r + wi*h2i;
	    *(buf+i4) = -h1i + wr*h2i + wi*h2r;
	}
	wr = ((temp = wr) * wpr) - (wi * wpi) + wr;
	wi = (wi * wpr) + (temp * wpi) + wi;
    }
    if ( forward )
	*(buf+1) = br;
    else
	[self cfft:forward];

    return self;
}

/* cfft replaces float array x containing NC complex values
   (2*NC float values alternating real, imagininary, etc.)
   by its Fourier transform if forward is true, or by its
   inverse Fourier transform if forward is false, using a
   recursive Fast Fourier transform method due to Danielson
   and Lanczos.  NC MUST be a power of 2. */

- cfft: (BOOL) forward
{
  float 	wr,wi,
		wpr,wpi,
		theta,
		scale;
  int 		mmax,
		ND,
		m,
		i,j,
		delta;

  ND = N2<<1;
  bitreverse( buf, ND );

  for ( mmax = 2; mmax < ND; mmax = delta ) {
    delta = mmax<<1;
    theta = twopi / ( forward? mmax : -mmax );
    wpr = -2. * pow( sin( 0.5*theta ), 2. );
    wpi = sin(theta);
    wr = 1.;
    wi = 0.;
    
    for ( m = 0; m < mmax; m += 2 ) {
      register float rtemp, itemp;
      
      for ( i = m; i < ND; i += delta ) {
	j = i + mmax;
	rtemp = ( wr * *(buf+j) ) - ( wi * *(buf+j+1) );
	itemp = ( wr * *(buf+j+1) ) + ( wi * *(buf+j) );
	*(buf+j) = *(buf+i) - rtemp;
	*(buf+j+1) = *(buf+i+1) - itemp;
	*(buf+i) += rtemp;
	*(buf+i+1) += itemp;
      }
      wr = ((rtemp = wr) * wpr) - (wi * wpi) + wr;
      wi = (wi * wpr) + (rtemp * wpi) + wi;
    }
  }

/* scale output */

  scale = forward ? 1./ND : 2.;
  
  { 
    register float *bi=buf, *be=buf + ND;
    
    while ( bi < be )
      *bi++ *= scale;
  }
  return self;
}

/* bitreverse places float array x containing N/2 complex values
   into bit-reversed order */

void bitreverse( float *buf,int N )
{
  float 	rtemp,itemp;
  int 		i,j,
		m;

  for ( i = j = 0; i < N; i += 2, j += m ) {

    if ( j > i ) {

      rtemp = *(buf+j);	/* complex exchange */
      itemp = *(buf+j+1);
      
      *(buf+j) = *(buf+i);
      *(buf+j+1) = *(buf+i+1);
      
      *(buf+i) = rtemp;
      *(buf+i+1) = itemp;
    }
    
    for ( m = N>>1; m >= 2 && j >= m; m >>= 1 )
      j -= m;
  }
}



/* oscillator bank resynthesizer for phase vocoder analyzer
   uses sum of N+1 cosinusoidal table lookup oscillators to 
   compute I (interpolation factor) samples of output O
   from N+1 amplitude and frequency value-pairs in C;
   frequencies are scaled by P */

- oscbank
{
  static int 	NP,
		tablesize = 8192,
		first = 1;
  static float 	Iinv,
		*lastamp,
		*lastfreq,
		*index,
		*table,
  	 	Pinc,
		ffac;
  int 		amp,
		freq,
		n,
		chan;

/* first pass: allocate memory to hold previous values
   of amplitude and frequency for each channel, the table
   index for each oscillator, and the table itself; also
   compute constants */

    if ( first ) {
  
      float twopioL = twopi / tablesize, tabscale;

	first = 0;
	lastamp = (float *) space( N2+1, sizeof(float) );
	lastfreq = (float *) space( N2+1, sizeof(float) );
	index = (float *) space( N2+1, sizeof(float) );
	table = (float *) space( tablesize+1, sizeof(float) );

	tabscale = ( Nw >= N2 ? N2 : 8 * N2 );

	for ( n = 0; n < tablesize+1; n++ )
	    *(table+n) = tabscale*cos( twopioL*n );

	Iinv = 1. / interp;
	Pinc = freqmlt * ( (float) tablesize ) / srate;
	ffac = freqmlt * pi / ((float) N2);
	if ( freqmlt > 1. )
	    NP = (int) ((float) N2 / freqmlt);
	else
	    NP = N2;
    }

/* for each channel, compute I samples using linear
   interpolation on the amplitude and frequency
   control values */

    for ( chan=0; chan < NP; chan++ ) {

      register float 	a,
			ainc,
			f,
			finc,
			address;

      freq = ( amp = ( chan << 1 ) ) + 1;
      if ( *(channel+amp) < synt ) /* skip the little ones */
	continue;
      *(channel+freq) *= Pinc;
      finc = ( *(channel+freq) - ( f = *(lastfreq+chan) ) ) *Iinv;
      
      ainc = ( *(channel+amp) - ( a = *(lastamp+chan) ) ) *Iinv;
      address = *(index+chan);
      
/* accumulate the I samples from each oscillator into
   output array O (initially assumed to be zero);
   f is frequency in Hz scaled by oscillator increment
   factor and pitch (Pinc); a is amplitude; */

      for ( n=0; n < interp; n++ ) {
	*(output+n) += a * *(table + ( (int) address ));
	address += f;
	
	while ( address >= (float) tablesize )
	  address -= (float) tablesize;
	
	while ( address < 0. )
	  address += (float) tablesize;
	
	a += ainc;
	f += finc;
      } 

/* save current values for next iteration */

      *(lastfreq+chan) = *(channel+freq);
      *(lastamp+chan) = *(channel+amp);
      *(index+chan) = address;
    }
    return self;
}


@end
